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This paper presents an asymptotic analysis of the Boltzmann equations (Riccati differential equa- 
tions) that describe the physics of thermal dark-matter-relic abundances. Two different asymptotic 
techniques are used, boundary-layer theory, which makes use of asymptotic matching, and the delta 
expansion, which is a powerful technique for solving nonlinear differential equations. Two different 
Boltzmann equations are considered. The first is derived from general relativistic considerations 
and the second arises in dilatonic string cosmology. The global asymptotic analysis presented here 
is used to find the long-time behavior of the solutions to these equations. In the first case the nature 
of the so-called freeze-out region and the post-freeze-out behavior is explored. In the second case 
the effect of the dilaton on cold dark-matter abundances is calculated and it is shown that there is a 
rsj . large-time power-law fall off of the dark-matter abundance. Corrections to the power-law behavior 

, are also calculated. 
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I. INTRODUCTION 



The thermal history of nonbaryonic dark-matter (DM) species is highly relevant to the shaping of the universe as 
we find it today. The existence of DM is based on evidence at many length scales. At the scale of galactic halos, for 
J-~H" example, DM explains the observed flatness of the rotation curves of spiral galaxies [lj. According to observations 
over the past twelve years, 23% of the energy of the universe consists of DM. This number has been obtained by 
best-fit analyses of astrophysical data to the Standard Cosmological Model, which is a Friedmann-Robertson- Walker 
cosmology involving cold DM as the dominant DM species. The modern data is based on observations of type-la 
supernovae [2| , the cosmic microwave background @, |j| , baryon oscillations @ , and weak-lensing data jf>] . It should 
be stressed that estimates of the DM abundance depend crucially on the theoretical model that is considered. 

Xj In the absence of dilaton effects from string theory, the evolution of the appropriately normalized number density 

qq ' Y(x) of a DM species X of mass mx is governed by the Boltzmann equation 



Y'(x) = -Xx- n - 2 [Y\x)-Y^{x)}, (1) 

which is a Riccati equation in the dimensionless independent variable x = mx /T, where T is the temperature. The 
parameter A is a dimensionless measure of the scattering of DM particles and is regarded as a large number A ^ 1. 
The integer n = 0, 1, 2, . . . comes from a partial- wave analysis of the scattering of DM particles; n — refers to 
S'-wave scattering. For bosonic remnants the function Y oq (x) is the distribution [7] 

x , 

_S: y^) = aj o ds eV ^_ i , ( 2) 

where A = 0.145g/.<7*, g is the degeneracy factor for the DM species, and g* counts the total number of massless 
degrees of freedom Q • 

As the universe cools and x increases, the nature of the solution Y(x) to ([lj changes rapidly in the vicinity of a 
value x = xj, the so-called freeze-out point, and as x — > oo the solution Y(x) approaches the constant Y^, called the 
relic abundance. Because a closed-form analytical solution to this Riccati equation is unavailable, an approximate 
heuristic approach is customarily used to treat this Riccati equation: One approximation is made for x < Xf and 
another is made for x > Xf. The solutions in the two regions are then patched at x = Xf. This approach gives an 
intuitive and reasonably accurate determination of Yoq and it is widely adopted [8| . 

However, this splitting into two regions is only a mathematical convenience and there is really no precise value 
Xf. Because the differential equation ([T]) is first order, its solution is completely determined by one initial condition, 
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namely Y(0). The usual method of splitting ([TJ into two approximate first-order equations, which are valid in each 
of two regions, leads to two conditions, an initial condition and a patching condition. We feel that this gives rise to 
an unsatisfactory mathematical discussion that is prevalent in the literature. The value of Xf, for example, becomes 
explicitly involved in the determination of Y^ when there is no reason for this. 

Equation (TT]) is valid in a general-relativistic framework. However, given the importance of understanding the 
current thermal-relic abundance of DM in theories beyond the standard model of particle physics, we also reexamine 
here the modifications of (jTJ) due to string cosmology [9]. String theory is widely accepted as a leading candidate 
for physics beyond the standard model, and it places a constraint on the types of time-dependent backgrounds in 
conformally invariant critical theories. As before, we are interested in eras in which the temperature T satisfies 
mx > T > T , where T is the current temperature of the universe. String cosmology leads to a rolling dilaton source 
in the Boltzmann equation [lOj that describes DM species. Including this source gives an additional linear term in 
the Boltzmann equation: 

Y'(x) = -Ax— 2 [Y 2 (x) - Y c 2 q (x)} + %Y(x)/x, (3) 

where <&o is a negative dimensionless constant of order 1 . 

The purpose of this paper is to study analytically the two Riccati equations ([T]) and ([3|). These equations do 
not have exact closed-form solutions. However, because A is a large parameter, one can attempt to find asymptotic 
approximations to the solutions. The most direct approach is to convert these Riccati equations into equations of 
Schrodinger type. When this transformation is applied to (fTJ, we obtain 

in \ n(n + 2) 2 _2n-4v2 / \ n I a\ 

V ( X > |^2 V \ X > ~ A X Y e q (X)V = 0. (4) 

Now, if we set n = 0, we obtain the standard time-independent Schrodinger equation in which 1/A plays the role of H. 
While it is possible to perform a local analysis of this equation for small x and for large a;, it is not easy to use WKB 
analysis to find a global asymptotic approximation because the equation is singular at x = and there is a turning 
point at x — oo. 

Thus, in this paper we will use two other powerful asymptotic methods from which we can extract global information. 
The first method is boundary-layer analysis. This asymptotic technique, which has been used to solve approximately 
the equations of fluid mechanics, gives very accurate results, and it has the physical advantage of treating freeze-out as 
a boundary-layer region, very much like the boundary between two fluids. The second technique, known as the delta 
expansion [111 ] , is particularly well-suited to study the transition from the equilibrium region to the large- a; behavior 
of the solutions without the necessity of finding approximations to the Boltzmann equation in different epochs. We 
will see that the presence of a source in ([3]) gives a solution for Y{x) in §5§, whose qualitative behavior is significantly 
different from the solution for Y(x) in JTJ). 

This paper is organized as follows. In Sec. [II] we summarize the derivation of the Boltzmann equations ([TJ and ([3]). 
In Sec. [HI] we a PPb/ boundary-layer analysis to study JT]) and ([3]). Next, in Sec. [IV] we describe the delta expansion 
and then use it to study the approximate behaviors of (Q]) and ([3]). Finally, in Sec. [V]we give some brief concluding 
remarks. 

II. DERIVATION OF THE BOLTZMANN EQUATIONS 

In this section we review the derivation of the two Boltzmann equations |T]) and ([3J . 

A. Derivation of ([l} 

In the hot early universe DM particles interact with themselves and with other particles. Particle species are assumed 
to react rapidly enough to maintain equilibrium. However, the universe expands and cools throughout its history. The 
timescale associated with this expansion is determined by the Hubble rate H. There is also a timescale T associated 
with the scattering cross-section (that is, an interaction rate per particle). The dynamics of DM particles depends on 
the ratio T/H. When T/H ^> 1, conditions for equilibrium hold and Y{x) follows the canonical distribution obtained 
from equilibrium statistical mechanics. However, for T/H <C 1 the DM particles are unable to maintain equilibrium. 
There is a crossover to freeze-out behavior in which Y(x) is asymptotically a constant. 

Let us consider a two-body scattering process in which particles of species 1 and 2 scatter reversibly into particles 
of species 3 and 4. The phase-space distribution function fi (r,p,t) for the species i gives the number of particles 



in an infinitesimal region of phase space around the position r and momentum p: fi (r, p, t) d 3 r d 3 p. The main bulk 
quantity of interest is the number density n^ (r , t) , which is given by [8j 



d 3 p 
(2nj~ 3 



n-i (r, t)=gi / 7—3 /; (r , p, t) , (5) 



where gi is the degeneracy factor for the ith DM species. The evolution of such a bulk quantity in the universe is 
given by the Liouville equation (in the absence of collisions, for simplicity) 



f-MI+f^+f 



^■ = mi=\i; + -§-V f +^-V f )f, = 0. (6) 



The standard Robertson- Walker metric for an isotropic and expanding flat universe is given by 

ds 2 = -dt 2 + a 2 (t) (dx 2 + dy 2 + dz 2 ) , (7) 

where a(t) is the scale factor [12J]. The covariant generalization of © is [8], 

where the Christoffcl symbol is given by 

^ vp = 9 \9av,p + 9ap,v 9vp,a) /*• 

For the metric in (|7J| , isotropy further implies that fi (p, t) — fi ( \p\ , t) . For the isotropic case ([5]) takes the form 



where _E = \J p 2 + m 2 . 

For a two-body collision process the Liouville equation ((SJ) no longer has a vanishing right side. This equation can 
then be used to describe the change in the number density of a given species. For species 1, for example, one gets Q 



a' 3 



d(nia 3 ) f d 3 pi f d 3 p 2 f d 3 p 3 f d 3 p 



dt J {2tt) 3 2E 1 J (2tt) 3 2E 2 J (2tt) 3 2E 3 J (2tt) 3 2E 4 

x (2^) 4 ,5 3 (pi +P2-P3- Pi) S (E 1 +E 2 -E 3 - E A ) \A\ 2 

x [f 3 U (1 ± fx) (1 ± h) - fxh (1 ± h) (1 ± U)] , (9) 

where the plus sign is used for a bosonic species and the minus sign is used for a fermionic species. The symbol A 
represents the scattering amplitude for the process 1 + 2 <-> 3 + 4 and it is a function of the pi. 

If the scattering process is sufficiently fast, fi can be parametrized by canonical Fermi-Dirac or Bose-Einstein 
distributions. For temperatures T <C E — /i the Bose-Einstein and Fermi-Dirac distributions both take the form 

f(E) ~ e^e-*/^ (10) 

which implies that quantum statistics are not important. Hence, the Pauli-blocking and Bose-enhancement are 
negligible (fi -c 1), and the third line of ([9]) simplifies: 

hf A (1 ± fl) (1 ± / 2 ) - /l/ 2 (1 ± h) (1 ± A) ~ e^ E - +E ^' T [ e (^+M4)/T _ e ( Ml +M 2 )/T 

where the relation E\ + E 2 = E 3 + E4 has been used. Also, combining ([5]) and (TTOl) . we get 

The equilibrium number density in the absence of a chemical potential is denoted by n\ . Thus, 

-3 d 1 3\ (0) (0), \) n 3 n 4 n x n 2 \ 

a _(n ia)= ni^^M —_-—_, (11) 



where the thermally averaged annihilation cross-section (av) is given by 

( m ,\ = X f d3pi f ^P 2 f d3p3 f d3pi r -(E 1+ E 2 )/T 

{ ' ~ n^n^J (2tt) 3 2eJ {2n) 3 2E 2 J (2w)*2E 3 J (2ir) 3 2E 4 

x (2tt) 4 S 3 (p! +P2-P3- Pi) S (E 1 +E 2 -E 3 - E A ) \A\ 2 . (12) 

We now make the standard assumption Q that the predominant interaction of the cold DM species X of mass n%x 
is XX -H- 11, where I is a light particle in equilibrium. As a consequence, in (|11[) we can replace n% and n 2 by nx, 
where nx is the number density of the species X. Also, we replace and n$ and 714 by n\ . The resulting equation is 



a 



d 



— (n x a 3 ) = (av) 



dt 



H 0) ) ; 



n\ 



(13) 



We now define x = m/T and note that dx/x — —dT/T = da/a because T scales as 1/a. Thus, ^f = i?x, where 
the Hubble rate if = ^log(a). Since the cosmological era for DM production is radiation dominated, a(t) oc \ft. 

This translates into H — H m /x 2 with H m = 1.67g* rn x /i~np\tinck- It is known theoretically [8| that av ex u 2n with 
n = for s-wave annihilation and n = 1 for p-wave annihilation. Since (v) oc \/T, we have the parametrization 
(av) = croa^" for x > 3 [|. 

Finally, we introduce the dependent variable Y = nx/T 3 ex nxa 3 . Similarly, we define Y cq = n x /T 3 . We then 
obtain the Boltzmann e quat ion in ([1]), where A = aom x / H m ex mpianck/"ix, and this explains why A is a large 
dimensionless parameter [13l |. 

B. Derivation of © 

String theory can be formulated in nonflat backgrounds, which is necessary when considering cosmology. Here, we 
consider the world-sheet sigma-model approach for dilaton-based cosmologies [9j. In superstring theory the bosonic 
part of the supermultiplet with lowest energy consists of the following massless states: the graviton gu n, the spinless 
dilaton $, and the antisymmetric spin-one tensor Bmn- For expanding universes $ provides consistent time-dependent 
backgrounds. In such backgrounds the string sigma model on the world sheet S is given by [14| 



S„ = 



d 2 a 
47ra' 



^j 1 ^g M N(X)d a X M dpX N + B M N(X)e a0 d a X M dpX N + a'^f$(X)RW/2\ , (14) 



where X M are target space-time coordinates with M,N — 0, 1, . . . , 9, a a are the world-sheet coordinates with 
a, = 0, 1, 7 Q/3 is the world-sheet metric, 7 = Idet (7 a/3 )|, R^ is the Ricci scalar associated with j a/3 , and a' is the 
string slope. Expanding around a conformal flat background with the action S*, we can write S a as 

S a = S* +h l f d 2 aV i7 (15) 

where h l denotes the background fields {<?mjv, Bmn, <&} and Vi are associated vertex operators [14|. 

Short-distance singularities of the quantum field theory on the world sheet lead to renormalized couplings {h l R j 
and to dependence on the renormalization-group scale u |15j . Usually, this results in nonvanishing j3 functions: 
ft 1 = dh l R /dlog(i. To restore conformal invariance, these (3 functions must vanish. This leads to equations of motion 
satisfied by the background fields. The usual procedure is to consider an effective target-space action in the string 
frame that reproduces the equations of motion: 



2a' 4 



/ d 10 x VG e~* \R + (V$) 2 + 2c/ 4 [/($) - H 2 /12~\ , (16) 



where H 2 = H^ va H^ va , H^ va = d^B va + d v B afJi + d a B^, and the potential U($) has been introduced. With the 
help of duality symmetries it is possible to find analytic solutions for the time dependence of the dilaton field [9| . 
From (fT5|) it can be shown [l(| that in three spatial dimensions the energy density p of the DM species X satisfies 

±+3H(p + p)-^(p-3 P )=0. (17) 



We then assume that the thermal DM species X behaves like dust (that is, p = 0) and that the energy density of the 
DM is given by the simple formula p = rrixnx- Next, in place of the on the right side of (|17p . we include a collision 
term, which is just the right side of (|T3l) : 



-n x + 3Hn x —n x = (av) 



(0) 

l x 



'■x 



(18) 



Assuming that matter sources are perfect fluids and requiring scale-factor duality symmetry, one can show [9J that 
up to an additive constant, $(£) = <&o loga(t) where $o = 0(1) and $o < 0. 

Finally, we make the assumption that the behavior of the DM species is dominated by radiation so that a(t) ex t 1 / 2 
0|. As in Subsec. [II A[ we introduce the variables Y(x) and x and obtain the Boltzmann equation ©. 



III. BOUNDARY-LAYER SOLUTION TO © AND © 



In this section we show how to perform a boundary-layer asymptotic analysis of |T]) and ([3]) . The advantage of this 
analysis is that it provides a global picture of the cosmological development from the initial time to the present as 
described by the Boltzmann equation, and not just the physics of the equilibrium epoch or of the post-equilibrium 
epoch alone. It also establishes a framework to describe in a clear and natural way the region of rapid transition 
between these two epochs. A satisfactory description of this crucial transition is lacking in earlier treatments in the 
literature because the earlier analysis used patching (joining together two solutions to a differential equation at an 
arbitrary and fictitious point, which produces an elbow in the solution) rather than asymptotic matching [16j . 

One may wonder why an asymptotic procedure as powerful as boundary-layer theory should be used to solve a 
first-order ordinary differential equation as simple as a Riccati equation. A general Riccati equation 

y'(x) = a(x)y 2 (x) + b(x)y(x) + c(x) 



can be recast as a linear second-order equation, 

'a'(x) 



w' (x) 



a(x) 



b(x) 



w'(x) + a(x)c(x)w(x) = 0, 



w' (x) 



where y(x) = r"\ r \ ■ Furthermore, (|19p can be recast as a Schrodinger equation 

v"{x) + {p"{x)/p(x) - [b(x) + a'{x)/a(x)] 2 /2 + a{x)c(x)} v{x) = 



(19) 



(20) 



by introducing v{x) = w(x)/p(x), where p'(x)/p(x) = [b(x) + a'(x)/a(x)]/2. 

The form (|2"0")) is often useful for asymptotic WKB analysis but the problem of freeze-out poses mathematical 
difficulties. If we apply these transformations to ((T|) for the case n = and use the leading asymptotic forms for 
y e q(^) in Appendices A and B, we obtain the Schrodinger equations 



where r\ = 2A£(3) for i<1, and 



v"(x) — A x T] v(x) = 



v"(x) - \ z A z x~ v e~ lx v(x) = 



(21) 
(22) 

(23) 

where v+and V- are constants. The approximate general solution to (|22l) can be found by using a standard application 
of WKB [I(|. [A detailed analysis of (|2"21 for large x is given in Appendix O] However, because (|2"2"|) has a turning 
point at x — 00 and because (|2ip has a singularity at x = 0, it is very difficult to construct a uniform asymptotic 
expansion that is valid for all x. We show below that boundary-layer theory overcomes these difficulties. 



for x 3> 1. The role of H in these equations is played by 1/A because A is treated as a large parameter. 
The exact general solution of (|2~TJ) is 



v{x)=x( K v + e Xrl/x +v-e- xrl/x ^j 



A. Boundary-layer analysis of (TjQ) 

Whenever the highest-derivative term in a differential equation is multiplied by a small parameter, one can attempt 
a boundary-layer analysis [161 ] . In such an analysis one identifies an outer region (or regions) in which the solution is 
slowly varying and an inner or boundary-layer region (or regions) in which the solution is rapidly varying. If these 
regions have an overlap, one tries to construct a global asymptotic approximation to the differential equation by 
performing an asymptotic match of the outer solutions to the inner solutions. 

In boundary-layer form the derivative term in (JT]) is multiplied by 1/A, which is regarded as small (1/A <C 1). Thus, 
we begin by looking for an outer solution; that is, a solution whose derivative is not large. To leading order such a 
solution in the outer region satisfies a distinguished limit (an asymptotic balance between two of the three terms in 
the differential equation) in which we neglect the derivative term as A — > oo: 

Y(x) ~ Y cq (x) (A -► oo). (24) 

Since Y(x) is well approximated by Y eq (x) in this region, we call this outer region the thermal- equilibrium region. 

To higher order, we seek a series expansion of this thermal-equilibrium outer solution as a formal power series in 
inverse powers of A: 

oo 
ythermal— equilibrium/ \ \ ' > — fcythermal— equilibrium/- \ (OK] 

fe=0 

Substituting this series into (Q} and collecting powers of 1/A yields the higher-order terms in the outer series. For 
example, to first order we get 

^thermal-equilibrium^ = _ 1 „+2 * ^ ^ cq (x)} . (26) 

We must now determine the extent of the thermal-equilibrium region. We know from Appendix A that for large x, 
i>1, the asymptotic behavior of Y oq (x) is given by 

Y cq (x) ~ Ae~ x x 3/2 (ar-»oo). (27) 

Thus, for large x in the outer region 



Y, 



thermal— equilibrium 



{x)~Ae- X X 3 ' 2 and y^mal -equilibrium^ „ i x n+2 (2g) 

Hence, the second term in the outer series is no longer small compared with the first term when 

a;~log(2A\)-(n+l/2)log(a;). (29) 

We will call the solution to this asymptotic relation the so-called freeze-out value Xf : 

xi ~ log(2AA) - (n + 1/2) log (a* ) . (30) 

Note that if we take A « 10 14 and A « 0.00145, we see that the outer asymptotic approximation ceases to be valid 
when x exceeds the approximate numerical value 

x { « 25. (31) 

Equation (|29j) defines the upper asymptotic limit of the thermal-equilibrium region. However, it is important to 
emphasize here that freeze- out does not occur at a point; Xf should not be viewed as a number but rather as a large 
range of values of x all satisfying the asymptotic relation (l29l) : 



x~xi (A->oo). (32) 

A second possible distinguished limit of (JlJ could in principle consist of an asymptotic balance between the left side 
and the second term on the right side. However, this distinguished limit is inconsistent and must be rejected because 
we are led to a contradiction: If we solve the resulting equation, we find that for large A the first term on the right 
side is in fact not negligible compared with the second term . 



A third distinguished limit of (JXJ) occurs when x is so large that the contribution of the equilibrium term Y£ (x) is 
negligible. In this case, the left side is asymptotic to the first term on the right side: 

Y'{x) ~ ~\x- n - 2 Y 2 {x) (x > 1). (33) 

In this second outer region, which we will call the post-freeze-out region, the solution yp^-fr 0020 - "*^) to (133")) is 

-impost— freeze - out /N _ ("1A\ 

(X) l/C-Xx-^/in+iy [ ' 

where C is a constant of integration to be determined. Note that this solution is consistent and valid when x ^> 1 
because Y oq (x) is exponentially small when x 3> 1. Note also that asx-} oo, Y post ~ h ' cczc ~ out (x) approaches the 
limiting value C. Thus, C represents the long-time limiting value of the relic abundance. 

The physical process of freeze-out can be recast in mathematical terms as a process that occurs in an inner region 
(or boundary layer), which we treat as a time interval that is comparatively short relative to the time intervals of 
the two outer regions, the thermal-equilibrium region and the post-freeze-out region. We begin the analysis of the 
freeze-out boundary layer by determining the size of this region. To do so, we introduce the inner variable X: 

x = Xf + kX. (35) 

We regard \X\ as a variable that may get large compared to 1, say as large as A max , but X is still small compared 
with A. Thus, since k is expected to be a small parameter roughly of order 1/A, the boundary layer is narrow because 
it extends roughly from X{ — nX max to X{ + KA" max . 

Making the change of variables (|3"5")) . from which we get 

— - IjL (w\ 

alx~ K dX' [6b) 

and treating kX as small compared with xi , we find that |T]) becomes 

-y'{X) = -Xx7 n - 2 \y 2 (X) - A 2 x\e~ 2xt ] , (37) 

K 

where y{X) = Y{x). A consistent dominant balance in this equation is achieved if we take 

k = x^ +2 /\, (38) 

and if we make this choice, we must neglect the second term on the right side because it is of order A -2 compared 
with the first term on the right side. This gives the simple inner differential equation 

y'(x) = -y 2 (x), (39) 

whose solution is 

y(X) ^ yijj. (40) 

where D is an integration constant. 

To complete the boundary-layer analysis, we must match the two outer solutions to this boundary-layer solution. 
In order to perform the asymptotic match, we re-express the outer solutions in terms of the inner variable X and 
then carry out an asymptotic approximation valid for small k to these asymptotic approximations. 

Let us look first at the outer solution in the post-freeze-out region: 

\^post— freeze— out / -v^\ /a -t \ 

1/C - A (xf + KX)-' 1 - 1 /(n + 1)' 
which simplifies to 



T^post — freeze— out / x^\ (A*)} 

X+ C ~ («+l)(x f )" + 1 



>s 



The coefficient of X in the denominator is 1, which agrees exactly with the coefficient of X in the inner solution (J40I) . 
Thus, we have achieved an asymptotic match, and the matching condition relates the constants C and D: 

D =h-i^hr- (43) 

Next, we match the boundary-layer solution in (|40[) to the outer solution (|25[) in the thermal-equilibrium region. 
To do so, we must re-express the outer solution in (|25[) in terms of the inner variable X. Although we are matching 
to just one term of the inner freeze-out solution, it is essential that we take the first two terms in the outer thermal- 
equilibrium series, and not just the first term, because we have shown that as we approach the freeze-out region, the 
first two terms in the outer solution become comparable in size. Thus, we include a factor of two in the asymptotic 
behavior 

^thermal— equilibrium/ \ ^ cy \ 3/2 — x 

~ 2A{ Xi + nXf /2 e- Xi e- KX 

- — ^- (44) 



X 



Because the coefficient of X in this behavior is 1, we obtain once again a perfect asymptotic match to the inner 
freeze-out solution in (14*01) . This allows us to determine the value of the constant D: 

D = \xj n - 2 . (45) 

Finally, combining this result with (|43| , we obtain the value of C: 



C = p^L, (46) 

A (n + 1 + Xf ) 

which is our result for the thermal-relic abundance. For Xf large compared with n + 1 this is in close agreement with 
the value (n + l)a;J l+1 /A given in Ref. 0. 

B. Boundary-layer analysis of ([3J 

The arguments given in Subsec. IIII Al apply to a modified version of (j3|). We modify (j3|) as follows. If we let 
ip = |$o|j then the substitution 

Z(x) = Y(x)x v 

reduces ((3]) to the simpler Riccati equation 

Z'(x) = -Xx- n - 2 [x- v Z 2 (x) - ^r c 2 q (a;)] . (47) 

The advantage of this equation over ([3]) is that there are only three rather than four terms, and thus it is easier to 
identify a dominant balance. 

We can now analyze (|47|) using the procedure adopted in the previous subsection. In the left outer region (the 
thermal-equilibrium region) we have 

^thermal-equilibrium^ _ ^-^+3/2 and ^thermal-equilibrium^ _ I^+n+2 _ ( 4 g) 

From this result we deduce that the freeze-out value Xf is given by 

xf - log(2^A) - (n + 1/2) log (a*) , 

which is identical to the result in (|30|) . This result shows that to leading order in 1/A the freeze-out temperature is 
independent of $o; that is, the location of the freeze-out region is only weakly affected by the presence of a dilaton. 
Next we discuss the right outer region (post-freeze region). The analog of ([34]) is 

ypost-freezc-out/ \ ( AQ\ 

[> i/c-\ x - n - i -*/{n + i + v y l ' 



where C is a constant of integration to be determined by asymptotic matching. As before, C describes the long-term- 
abundance behavior. However, when x is large compared with the freeze-out temperature (x 3> Xf), Y(x) does not 
approach a constant. Rather, 

Y(x) ~ x~ v Z{x) ~ x~ v C (x -> oo). (50) 

In the freeze-out boundary-layer region we again make the change of variable in (|35p , 

x — Xf + nX, 

where the inner variable X may become large compared to 1, but it is still small compared with A. Thus, since k is 
expected to be a small parameter of order 1/A, the boundary layer is narrow as before. A consistent dominant-balance 
gives the value 

« = x* +2+v /\. (51) 

The inner differential equation then has the form 

Z\X) = -Z\X\ (52) 

where Z(X) = Z(x). The solution to ([52]) is 

Z(X) = _L_, (53) 

where D is an integration constant. This is the analog of (p7)l . 

An asymptotic match of the right outer solution to the boundary-layer solution produces the relation between the 
constants C and D, 

o (n + 1 + (p)x { 

which is the analog of (|43| . Finally, by matching the left outer solution to the boundary-layer solution, we obtain the 
value of C: 

C='" + 1 r W "T (55, 

A (n + 1 + <p + xi ) 

In conclusion, we find that, due to the presence of a dilation, the thermal-relic abundance in (J50I) remains time 
dependent; it vanishes as x — > oo and does not approach a constant. Note also that if we eliminate the effect of the 
dilaton by allowing $o to approach 0, the results in ((511)) and ([53)1 smoothly reduce to that in Subsec lHI A[) . 

IV. APPLICATION OF THE DELTA EXPANSION TO (© AND © 

In this section we show how to apply the delta expansion to (QJ and <j3j) . We begin with a brief summary of the 
delta-expansion technique. 

A. Summary of the delta expansion 

The delta expansion is an unconventional perturbative technique for solving nonlinear problems. It was first 
introduced to treat nonlinear aspects of quantum field theory |17| . To prepare for applying it to the Boltzmann 
equations (P) and §5$ , in this subsection we give a brief review of the delta expansion. 

The theme of the delta expansion is to introduce a parameter 5 as a measure of the nonlinearity of a problem; that 
is, the departure of the problem from a corresponding linear problem. We then treat S as small (6 <C 1), and solve 
the problem perturbatively by expanding about the linear problem obtained by setting S = 0. The basic ideas of the 
delta expansion are explained in Ref. [111 ]. 

To illustrate the delta expansion, we consider the Thomas-Fermi nonlinear boundary-value problem 

y"(x) = [y(x)] 3 / 2 /^, 2/(0) = 1, y(+^) - 0. (56) 
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This problem is extremely difficult and no closed-form analytical solution is known. We introduce the parameter S in 
the exponent of the nonlinear term of the differential equation and consider the one-parameter family of problems 

y"(x) = y(x)[y(x)/x] s , y(0) = 1, y(+oo) = 0, (57) 

where we treat 8 as a small perturbation parameter. The solution to the unperturbed (8 — 0) linear problem is 
yo(x) = e~ x , and we use y^ix) as the first term in the delta expansion of the solution to the nonlinear problem (1571) : 



i{x) = Y J 5 k y k {x). (58) 



fc=0 

Finally, we recover the solution to the original Thomas-Fermi problem by setting 8 = 1/2. Typically, only very few 
terms are needed in the delta expansion to recover accurate numerical results. Furthermore, the accuracy of the delta 
expansion can by accelerated by using Pade techniques to sum the delta expansion. In the case of the Thomas- Fermi 
problem a (2, 1)-Pade approximant has a numerical error of about 1%. 
As a second example, consider the quintic polynomial equation 

x 5 + x - 1 = 0, 

which cannot be solved by quadrature. The real root of this equation is x — 0.75487767 . . .. Introducing the 
perturbation parameter 8, we obtain the equation 

x 1+s +x = 1. 

We then seek a perturbation series of the form 

x(S) = c + ad + c 2 8 2 + c 3 S 3 + ... (59) 

whose first term is cq = 1/2. The radius of convergence of the delta series (J59I) is 1, and therefore it diverges at 8 = 4. 
However, a (3,3)-Pade approximant has a numerical error of 0.05% and a (6, 6)-Pade approximant has a numerical 
error of 0.00015%. 

B. Delta expansion for ([l]) 

To apply the delta expansion to ([1]), we insert the parameter 8 in such a way that when 8 = 1 we recover ([Tj): 

Y'(x) = -Xx- n - 2 (Y - Y cq ) (Y + Y cq ) s . (60) 

There are, of course, many ways to insert the parameter 8, but the advantage of (J60I) is that the solution to the 
unperturbed linear problem obtained by setting 8 — is qualitatively similar to the solution to ([T]), which we have 
already investigated in Sec. IIIH In particular, when 8 = 0, Y(x) behaves like Y cq (x) for small x, undergoes a transition 
as x increases, and then approaches a constant as x —¥ oo. 

Following the usual delta-expansion procedure, we represent Y(x) as a series in powers of 8, 

oo 
k=Q 

and then substitute this series into (|60p . Comparing powers of 8, we obtain a sequence of inhomogeneous differential 
equations for y^: 

y' k (x) + Xx- n - 2 y k (x) = h k {x) (k = 0, 1, 2, . ..), (61) 

where 

h (x) = Xx- n - 2 Y cq (x), 

hx{x) = Ax-"- 2 [y oq (a;)-yo(a;)]log[ycq(^)+2;o(a;)], 



h 2 (x) = \x- n - 2 yi {x) ^\ *" + Lqv ', " uw log^ [Y oq {x) + y (x)] - yi (x) log [y oq (a;) + y Q (x)] , (62) 
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and so on. 

The solution to (|5T|) . which is obtained by using the integrating factor exp [— \x~ n ~ l / (n + 1)] , has the quadrature 
form 



Vk{x) 



- „Ax-"- 1 /(n+l) 



ds e- Xx 



~ 1/{n+1) h k (s). 



(63) 



Because (|61|) is a first-order equation, its solution contains one arbitrary constant for each k and this constant is 
determined by the requirement that 2/fc(0) be finite. This requirement fixes the lower endpoint of integration to be 
for all k. Note that if we evaluate the integral in (j6"3")l . we obtain the results yo(0) = ^cq(O) = 2A((3) (see Appendix 
B), yi (0) = j/2(0) = . . . = 0. As a; increases, yo(x) remains close to Y oq (x) until x is of order A. 

We can now express the freeze-out value Y(po) as a series in powers of S and then evaluate this series at S = 1. 
Here, we just calculate the first term in ths series: 



y (x) = Ae^-VCn+D 



dss- n - 2 e- Xs 



" /(n+1) r eq ( s ). 



(64) 



Let us evaluate this integral assuming that the parameter A is large. Since the integrand is exponentially small for 
small s, we may assume that the only contribution to the integral comes from the region s> 1, and in this region we 
may replace 5^ q (s) by its asymptotic behavior As 3 ^ 2 e~ s (see Appendix A). We thus obtain 



y (oo) ~A\ d8»- n-1/ V w (A -> 00), 



(65) 



where 



4>{s) 



To evaluate (|65|) we use Laplace's method with a moving maximum [161 ]. We note that the maximum of </>(s), which 
occurs when 4>'{s) = 0, is at sq = \ 1 ^ n+2 \ Hence, we introduce the rescaled variable t: 



= u l/(„+2)_ 



This gives the integral 



where 



,'/o 



(00) ~ AA 5 /( 2 "+ 4 ) / dtr"- 1 /V 1/(B+!0 *tt (A ->oo), 
Jo 



(66) 



6{t) 



■t- 



1 



-t 



-n— 1 



1 



The maximum of 9(t) occurs at t = 1, and near this point we have the quadratic approximation 

n+2 n+2 



6{t) 



n+l 2 

Thus, evaluating the Gaussian integral, we obtain the result 

AV2^ 



-(t-iy 



2/0(00) 



\M + 2 



A 2/(n+2) 



exp 



1 



:a 1 /(«+2) 



which reduces to 



2/o(oo) ~ AX^pHe 



-1sf\ 



(67) 



(68) 



when n = 0. Thus, the delta expansion predicts that at x = 00 the freeze-out value of Y(x) is exponentially small. 

We see from this calculation that the delta expansion gives a simple and qualitatively accurate picture of the 
solution to the Boltzmann equation (Q}. However, the prediction in (|67| of the relic abundance YJx, is clearly too 
small and, of course, this is because we have only kept the leading-order term in the delta expansion. We will see in 
the next subsection that if we retain higher powers of 5, the qualitative features of the solution do not change but the 
quantitative prediction for the long-time behavior of Y(x) is improved. 
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C. Delta expansion for ([3]) 

The delta expansion treatment of ([3]) parallels that for ([l). We insert the parameter 5 into ([3]) as follows: 

A 



Y '( x ) = -~-^[Y(x)-Y cq (x)}[Y(x)+Y cq (x)} 5 - Z Y (x% (69) 

where <f> = |$ |- The analog of (f6"T|) is then 

J^*) + (^ + ^) »*(*) = &*(*) (fc = 0, 1,2, ...)• (70) 

The solution to (|70]). which is obtained by using the integrating factor x^ exp [~Aa; _ ™ _1 /(n + 1)] , has the quadrature 
form 

y k (x) = x-' t 'exp[Xx- n - 1 /(n + l)] f dss^exp [-Xs- n ~ 1 /(n + 1)] h k (s). (71) 

Jo 

Using the modified Laplace method again, we obtain for x — > oo and large A the asymptotic approximation 

y (x) ~ x~*B{\), (72) 

where the constant B(X) is given by 



B(X) = AA^+^W^exp [-A 1 /'^ 2 ) (n + 2)/(n + 1) 



This shows that the dilatonic correction to the Boltzmann equation gives a significant qualitative change in the freeze- 
out behavior of DM. The magnitude of the DM abundance is era dependent because its leading behavior for large x 
is an algebraic decay of the form x~^. The delta expansion is qualitatively in agreement with boundary layer theory. 
The result in (|72|) is the analog of (|67|) . and again we see that while the delta expansion in leading-order gives a 
good qualitative description of the solution to the Boltzmann equation, the quantitative prediction for the coefficient 
B(X) of x~^ in the large-x behavior is much too small. Thus, we extend the result in (|72|) to first order in 5. The 
calculation is a straightforward generalization of the zeroth-order calculation and the result is 



2/o (x) + S yi (x) ~ x-*B(\) 1 - <Jlog[B(A)] + 5 



n+1 



, ■ A 

7 + lo; 



n + 1 



(73) 



where 7 = 0.5772 ... is Euler's constant. 

For large A, we can ignore all but the log[i?(A)] term, and we obtain a rough asymptotic behavior, which is a 
simplified version of ([75)1 : 

yo(x) + S yi (x) - x~+B(\) {1 -5\og[B(X)}} . (74) 

Not surprisingly, the second-order contribution contains a logarithm squared: 

yo(x) + S Vl (x) - x^B(X) |l - *log[B(A)] + \f log 2 [5(A)] | . (75) 

In general, the dominant contribution to the coefficient of 5 k in the delta expansion is (— l) fe log [B{X)]/k\. Thus, 
if we sum the approximate delta series to all orders in 5 and set 5 = 1, the multiplicative coefficient B{X), which is 
numerically incorrect because it is much too small, is exactly canceled. This explains the mechanism by which the 
delta expansion and the matched asymptotic analysis become compatible. 

V. BRIEF CONCLUDING REMARKS 

We have applied two powerful perturbative techniques, boundary-layer theory and the delta expansion, to find 
globally accurate solutions to two different Boltzmann equations that describe dark-matter abundances in the early 
universe. The first Boltzmann equation is based on the standard model of particle physics and general relativity; 
the second includes additional effects due to dilatonic contributions that arise in string theory. The boundary-layer 
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solution consists of contributions from three distinct eras, a thermal-equilibrium epoch, a freeze-out region, and a 
nonequilibrium relic-abundance epoch, and the global solution is obtained by the use of asymptotic matching. The 
delta-expansion solution does not require the use of asymptotic matching and gives a good qualitative picture of the 
behavior in these three epochs, but the results to low orders in <5 are not as accurate for long times. 

We have shown that when dilatonic effects are not included, the dark-matter-relic abundance approaches a constant 
for long times, but when dilatonic effects are included, the relic abundance has a power-law decay determined by the 
dilaton coupling. 

Acknowledgments 

We are grateful to N. E. Mavromatos for enlightening discussions. CMB thanks the U.K. Leverhulme Foundation 
and the U.S. Department of Energy and SS thanks the U.K. Science and Technology Facilities Council for financial 
support. 

Appendix A: Large-x behavior of Y cq (x) 

In this Appendix we derive the large- a; asymptotic behavior of the equilibrium distribution Y cq in ([2]), whose integral 
representation is given by 



Ye « {x)=A Lo dS e^-l (A1) 

When x ^> 1, we can neglect —1 in denominator of the integrand to all orders in powers of 1/x and write 

/•OO 

Y cq (x)~A dss 2 e^ /sA+xi (ar-)-oo). (A2) 



The scaling s — xt followed by the change of variables u = \/t 2 + 1 then gives the integral representation 



Y cq (x) ~ Ax 3 / due- xu uVu 2 -l (x -> oo). (A3) 



Watson's lemma |16| applies directly to the integral (|A3|) . The procedure is first to expand uy/u 2 + 1 as a series in 
powers of u — 1, 



iV« 2 -l = $>«(«- l) n+1/2 , (A4) 



n=0 

where 



a n =^-l,2)^ + ^ n -'/ 2 \ (A5) 



/2tt * n\ 

and then to interchange orders of summation and integration. Integrating term by term gives the asymptotic series 

V27T „ (-2x) n n\ 

Thus, the series begins 

Y cq (x) ~ Ae^x 3 / 2 ^ (l + g + . . .) (x -► oo). (A7) 
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Appendix B: Small-a; behavior of Y cq (x) 

In this appendix we show how to find the small- x asymptotic behavior of the integral 



Y eq (x) = A da—== . (Bl) 

=o e ^ s + x - 1 



We begin by substituting t = y/s 2 + x 2 . This gives 



< x ,wt 2 -x 2 t r°° „ t 2 

't=x 

where 



Y cq {x)=Al dt v t _ =A dt^—-(^l-x 2 /t 2 -l + l\ = A + B + C + V, (B2) 



r°° t 2 

A = r oq (0) = A / dt^— = 2A((3), 
Jt=o e — 1 



r x t 2 
B = -A dt— 

Jt=o e 1 

r°° +2 , . 

c = ^ i ^- F - T (yi-^ W _ 1 ), 

2? = A / dt 

e 



i— a: 



^-(Vl-^A 2 -!)- (B3) 



We now evaluate each of the integrals B, C, and D, in turn. 

To evaluate B we expand t/(e t — 1) in a Taylor series, which converges if |i| < 27r, and integrate term by term: 



f' x °° D °° p 

B = -A dtty ^t n = -AY - \ x n+2 , (B4) 

Jo f^ n n\ f- n + 2k! 



where B n is the nth Bernoulli number (B Q = 1, B x = —1/2, B2 = 1/6, -B4 = —1/30, . . ., B 2n +i = for n > 1). So, 



To evaluate C and I? we use the expansion 



^.,. >j;!fi^. (B6) 

2v7r Z -^ n! 

v n— 1 



Thus, C becomes 

C 2^A e*-l^ n! 2^ (n + 1)! A *e*-l' (B7J 

Hence, 

C = c 2 x 2 + c 4 x 4 + O (x 6 ) , (B8) 

where 

1 f°° dt 1, „ 1M , l [ x dt 

02 = -^/ ^r-2 l0g(l - l/e) and c ^sl w^ry (B9) 

The interesting contribution comes from D. We express D as the double sum 

p = - * y ^ y r( / n+l {, 2) x 2 -+ 2 /' dtr- 1 - 2 ". (bio) 

v m= n=0 V ; Ji - 1 
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Depending on the values of to and n in the sum, we get different kinds of terms. For example, logarithm terms appear 
when (and only when) to = 2n. Thus, all the logarithm terms appear in the series 

V,, _ - ^L £ ^ ff ^» log* - ^ log* - 1,' log, + . . . . (Bll) 

Terms of order a; 2 arise from the upper endpoint of integration in (|B10p when n = and for all to > 1 (but not 
to = because this gives rise to a log term, and we have already included this contribution) and they arise from the 
lower endpoint of integration when m = for all n > 1 (but not n = 0). The upper endpoint gives 

Supper, 2 = ~\x 2 J dt (-^-j - ij = -1 l0g(l - l/e)x 2 . (B12) 

The lower endpoint gives 



IV„, , ,> = x 2 / L | [Jl-\t* + 1? -U=- A x 2 - 1 - log(2)x 2 . (B13) 



i=0 i 3 V V 2 A 2 



Thus, the result for the expansion of I Qq (x) in (|B1[) to order x is 

!(«) ~ 2C(3) 



^log(z/2)-i 



.x 2 + . . . . (B14) 



We have pursued this calculation to higher order in powers of x, and we find that in (|B14[) the coefficients of a; 3 
and a: 5 are 0, the coefficient of x A is 

7 Ct}l _ 1 _ lof(2) logjz) = _ 00297+001041 (x) (B15) 

96 8 128 96 96 6W v y 

and the coefficient of x 6 is 

7 C ; (-l) 979 7T log(x) 

192 16 268800 2880 11520 

where 7 = 0.57721 is Euler's constant. 



= -0.0048 + 0.0000868 log(x), (B16) 
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